#include<stdio.h>
#include<math.h>
#include<string.h>
#include<stdlib.h>
#include<ctype.h>
#include<dirent.h>
#include<glob.h>
#include <unistd.h>


namespace qi = boost::spirit::qi;

/*
rabid:/scratch/xiang/ima_new/magicpool/bcarab_1st/fqs/*.tmp2
rabid:/scratch/xiang/ima_new/magicpool/bcarab_2nd/fqs/*.tmp2
rabid:/scratch/xiang/ima_new/magicpool/bcarab_3rd/recon/*.tmp2

reconstruction -a /data/www/18genomes/variants.tables -m ./transposons.gff -c /Net/rabid/data/scratch/xiang/ima_new/magicpool/bcarab_1st/fqs/

*/



//!_______________________________________________________________________________________
//! 
//!		read_allele_data
//!			
//! 
//!		\author lg (05/04/2011)
//!     
//!     \param filename 
//!     \param alleles_per_contig 
//!_______________________________________________________________________________________
void read_allele_data(const string& file_name, vector<t_contig_alleles>& alleles_per_contig)
{
	opt.both() << "Reading from " + file_name + "\n";
	io::stream<io::mapped_file_source> input_stream;
	input_stream.open(file_name.c_str());
	if (!input_stream)
		throw std::logic_error("Can't open " + file_name);



	// progress dots
	t_progress_indicator progress(opt.verbose(), guess_rows / 100 + 1, guess_rows);
	progress.use_timer();


	//	
	//	Add data for new contig
	//
	alleles_per_contig.push_back(t_contig_alleles());
	t_contig_alleles& ad = alleles_per_contig.back()


	//
	//	Parse strain names from header line
	//
	string first_line_str;
	getline(input_stream, first_line_str)
	qi::phrase_parse(beg, end, +qi::graph , qi::space, ad.strain_names);
	if (ad.strain_names.size() <= 5)
		throw std::logic_error("No strain names in header of allele file [" + file_name + "]");
	if (ad.strain_names[0] != "chr")
		throw std::logic_error("Ill formed header line does not begin with chr: [" + file_name + "]");
	
	// get rid of first five fields which are before the strain_names
	ad.strain_names.erase(ad.strain_names.begin(), ad.strain_names.begin() + 5)
	ad.cnt_strains = strain_names.size();




	// 
	// ignore cases which are non-diallelic
	// 
	if ( nalleles != NULL && sscanf( nalleles, "%d", &Nalleles ) == 1 && Nalleles==2 )
	{
		sscanf( pse, "%d", &Pse );
		sscanf( bp, "%lf", &x );
		Bp = (int)x;
		sscanf( maf, "%d", &Maf );


		// save contig name 
		if ( n==1 )	ad->chr = strdup(chr);
		ad->all_alleles[n] = (char**)calloc(ncol,sizeof(char*));
		ok = 1;
		for (k=0;k<ncol;k++)
		{
			allele = strtok(NULL," \t");
			//
			//	if non SNP or numeric allele, OK = false 
			// 
			if ( (strlen(allele) > 1) || isdigit(allele[0]) ) ok *= 0;


			//
			//	copy first letter of each allele to buf 
			// 
			buf[k] = allele[0];
			//
			//	all_alleles = all saved allele 
			// 
			ad->all_alleles[n][k] = strdup(allele);
		}
		// 
		// 
		// calculate strain distribution pattern
		//  assigning indices to successively encountered new allels
		// 
		ad->sdp[n] = (char*)calloc(ncol,sizeof(char));
		ad->sdp[n][0] = 1;
		sdp = 1;
		for (k=1;k<ncol;k++)
		{
			int j;
			ad->sdp[n][k] = -1;
			for (j=0;j<k;j++)
			{
				if ( !strcmp(ad->all_alleles[n][k],ad->all_alleles[n][j]) )
				{
					ad->sdp[n][k] = ad->sdp[n][j];
					break;
				}
			}
			if ( ad->sdp[n][k] < 0 ) ad->sdp[n][k] = ++sdp;
		}


		//
		//	if any non SNP or numeric allele, OK = false 
		// 
		ad->ok[n] = ok;

		// 
		// save	SNPs to allele.
		// 
		//  	Note that first letter of non SNPs will also be saved but we
		//  	are guided by OK flag
		// 
		if ( ! ad->diallellic_snps[n] ) ad->diallellic_snps[n] = (char*)calloc(ncol+1, sizeof(char));
		for (k=0;k<ncol;k++)
		{
			ad->diallellic_snps[n][k] = buf[k];
		}

		// 
		// position along contig
		// 
		ad->loci[n] = Bp;
		n++;
	}
	nrow++;
}
}


	using qi::graph;
	using qi::int_;
	using qi::unsigned_;
	using qi::float_;
	auto grammar = +graph >> int_ >> float_ >> unsigned_ >> int_ >> +graph % qi::space;

	while (getline(input_stream, line_str))
	{
		++progress;
		cnt_rows += 1;

		// 
		// 	Parse whitespace separated elements
		// 		Name unused fields just in case...
		// 
		string 			contig;
		int				pseudo_pos;
		float			fractional_pos;
		unsigned		cnt_alleles;
		int				maf;
		vector<string>	alleles;
		auto beg = line_str.begin();
		auto end = line_str.end();
		//    chr   pse.bp   bp     nalleles   maf   SPECIES_1   SPECIES_2   SPECIES_3    
		//    1     11       11     2          18    T           T           T            
		//    1     22       21.5   3          16    0           0           0            
		qi::phrase_parse(beg, end, grammar , qi::space, fields, alleles);

		// check right number of alleles
		if (alleles.size() != ad.cnt_strains)
			throw std::logic_error(
									(format("Expecting %1$d alleles on every row. "
										   "Row %2$d has %3$d values.") 
										% ad.cnt_strains
										% cnt_rows
										% alleles.size()).str());
		int				pos = unsigned(fractional_pos);


		// 
		// save contig name
		//
		if (!ad.contig.length())
			ad.contig = contig;


	}

	if (!tmp_data_matrix_storage.size())
		throw std::logic_error(file_name + " is empty.");

	//
	//  copy 1 dimension vector into matrix
	// 
	data_matrix.resize(extents[cnt_rows][cnt_columns]);
	data_matrix.assign(tmp_data_matrix_storage.begin(), tmp_data_matrix_storage.end());
} 
}




typedef struct
{
	char *chr;
	int cnt_diallelic_snp_sites;
	int *loci;
	char *ok;
	int cnt_strains;
	int maximum_locus;
	char **strain_names;
	char **diallellic_snps;
	char ***all_alleles;
	int *map_pos_to_snp_index;
	char **sdp;
} ALLELE_DATA;

typedef struct
{
	int chr, from, to, id;
	double error;
} SEGMENTS;


//Chr1	11	T	?
//Chr1	18	T	?
//Chr1	73	C	?
typedef struct
{
	char **name;
	// size of  buffer
	int len;
	// position
	int *bp;
	int *chr;
	// 				SNP alleles
	char *allele;
	// pointer array of alleles in founder strains corresponding to this contig position
	char **founder;

	// list of founder strain names
	int Nfounders;
	char **founder_name;

	// contig bounds
	int *chr_start;
	int *chr_end;
	SEGMENTS **mosaic;
	int *nseg;
} SEQ_CALL;





ALLELE_DATA *ReadAlleleData( char *filename );
ALLELE_DATA **ReadGenomeAlleleData( char *dir, char *maskfile, int nchr );
SEQ_CALL **ReadSeqCalls( char *dir, char *suffix, ALLELE_DATA **add, int nchr, int *nseq );
SEQ_CALL *ReadSeqCall( char *filename, ALLELE_DATA **add, int nchr );
SEGMENTS *ReadGFF( char *filename, int *nsegments );
int reconstruct_mosaic_diploid( SEQ_CALL *sq, int penalty, int chr );
int reconstruct_mosaic( SEQ_CALL *sq, int penalty, int chr );


int main( int argc, char **argv ) {

	ALLELE_DATA **add;
	int index;
	int c;
	int diploid=0;
	int penalty=10;
	char *alleledir = NULL;
	char *calldir = NULL;
	char *maskfile = NULL;
	char *Penalty = NULL;
	opterr = 0;

	while ((c = getopt (argc, argv, "da:m:c:p:")) != -1)
		switch (c)
		{
		case 'd':
			diploid=1;
			break;
		case 'a':
			alleledir = optarg;
			break;
		case 'c':
			calldir = optarg;
			break;
		case 'm':
			maskfile = optarg;
			break;
		case 'p':
			Penalty = optarg;
			break;
		}
	if ( Penalty )
		sscanf( Penalty, "%d" , &penalty);

	if ( alleledir && calldir )
	{
		int s, chr, nchr=5, nseq=0;
		ALLELE_DATA ** add = ReadGenomeAlleleData( alleledir, maskfile, nchr ) ;
		SEQ_CALL **call = ReadSeqCalls( calldir , "tmp2", add, nchr, &nseq );
		for ( s=0;s<nseq;s++)
		{
			for (chr=1;chr<=nchr;chr++)
			{
				if ( diploid )
					reconstruct_mosaic_diploid( call[s], penalty, chr );
				else
					reconstruct_mosaic( call[s], penalty, chr );
			}
		}
	}
}


//!_______________________________________________________________________________________
//!     
//!		ReadAlleleData
//! 
//! 		allocates memory for allele and returns pointer
//!     
//!		\author lg (4/4/2011)
//!     
//!     \param filename 
//!     
//!     \return ALLELE_DATA* 
//!_______________________________________________________________________________________
ALLELE_DATA *ReadAlleleData( char *filename ) {

	FILE *fp = fopen( filename, "r");
	if ( fp == NULL )
	{
		fprintf(stderr, "ERROR Could not open file %s\n", filename );
		exit(1);
	}
	else
	{
		ALLELE_DATA *ad = (ALLELE_DATA*)calloc(1,sizeof(ALLELE_DATA));
		// read buffer
		char *buffer = (char*)calloc(1000000,sizeof(char));
		// genotypes
		char **strain_names = (char**)calloc(10000,sizeof(char*));
		int n=0, nrow=0, maximum_locus=0, k;
		int ncol=0;
		int Nbuf = 1000000;


		// n is the number of row with saved, allelic data 
		// nrow = physical number of rows
		while ( fgets( buffer, 1000000, fp ) )
		{
			/*      printf("%d %s\n", n, buffer); */
			buffer[strlen(buffer)-1]=0;
			if ( n==0 )
			{
				// 
				// read header row
				// 
				n++;
				char *chr, *pse, *bp, *nalleles, *maf, *allele;
				int Pse, Bp=0, Nalleles=0, Maf=0;
				chr = strtok( buffer, " \t");
				pse = strtok( NULL, " \t");
				bp = strtok( NULL, " \t");
				nalleles = strtok( NULL, " \t");
				maf = strtok( NULL, " \t");
				while ( allele = strtok( NULL, " \t"))
				{
					strain_names[ncol++] = strdup(allele);
				}
				// set up memory for all strain (column) names
				ad->strain_names = strain_names;
				// first row, read number of strains (columns)
				ad->cnt_strains = ncol;
				ad->cnt_diallelic_snp_sites = 0;
				ad->loci = (int*)calloc(Nbuf,sizeof(int));
				ad->ok = (char*)calloc(Nbuf,sizeof(char));
				ad->allele = (char**)calloc(Nbuf,sizeof(char*));
				ad->all_alleles = (char***)calloc(Nbuf,sizeof(char**));
				ad->sdp = (char**)calloc(Nbuf,sizeof(char*));
			}
			else
			{
				//    chr   pse.bp   bp     nalleles   maf   SPECIES_1   SPECIES_2   SPECIES_3    
				//    1     11       11     2          18    T           T           T            
				//    1     22       21.5   3          16    0           0           0            
				//    1     38       29.5   9          10    CTCTGA      CTC         CTCTGA       
				char *chr, *pse, *bp, *nalleles, *maf, *allele;
				int Pse, Bp=0, Nalleles=0, Maf=0, k, ok;
				char *buf = (char*)calloc(ncol+1,sizeof(char));
				double x;
				int sdp;

				chr = strtok( buffer, " \t");
				pse = strtok( NULL, " \t");
				bp = strtok( NULL, " \t");
				nalleles = strtok( NULL, " \t");
				maf = strtok( NULL, " \t");

				// 
				// ignore cases which are non-diallelic
				// 
				if ( nalleles != NULL && sscanf( nalleles, "%d", &Nalleles ) == 1 && Nalleles==2 )
				{
					sscanf( pse, "%d", &Pse );
					sscanf( bp, "%lf", &x );
					Bp = (int)x;
					sscanf( maf, "%d", &Maf );


					// 
					// Reallocate memory
					// in exponential way for efficiency
					// 
					if ( n >= Nbuf )
					{
						Nbuf *= 1.5;
						ad->loci = (int*)realloc(ad->loci,sizeof(int)*Nbuf);
						ad->ok = (char*)realloc(ad->ok,sizeof(char)*Nbuf);
						ad->allele = (char**)realloc(ad->allele,sizeof(char*)*Nbuf);
						ad->all_alleles = (char***)realloc(ad->all_alleles,Nbuf*sizeof(char**));
						ad->sdp = (char**)realloc(ad->sdp,Nbuf*sizeof(char*));
					}

					// save contig name 
					if ( n==1 )	ad->chr = strdup(chr);
					ad->all_alleles[n] = (char**)calloc(ncol,sizeof(char*));
					ok = 1;
					for (k=0;k<ncol;k++)
					{
						allele = strtok(NULL," \t");
						//
						//	if non SNP or numeric allele, OK = false 
						// 
						if ( (strlen(allele) > 1) || isdigit(allele[0]) ) ok *= 0;


						//
						//	copy first letter of each allele to buf 
						// 
						buf[k] = allele[0];
						//
						//	all_alleles = all saved allele 
						// 
						ad->all_alleles[n][k] = strdup(allele);
					}
					// 
					// 
					// calculate strain distribution pattern
					//  assigning indices to successively encountered new allels
					// 
					ad->sdp[n] = (char*)calloc(ncol,sizeof(char));
					ad->sdp[n][0] = 1;
					sdp = 1;
					for (k=1;k<ncol;k++)
					{
						int j;
						ad->sdp[n][k] = -1;
						for (j=0;j<k;j++)
						{
							if ( !strcmp(ad->all_alleles[n][k],ad->all_alleles[n][j]) )
							{
								ad->sdp[n][k] = ad->sdp[n][j];
								break;
							}
						}
						if ( ad->sdp[n][k] < 0 ) ad->sdp[n][k] = ++sdp;
					}


					//
					//	if any non SNP or numeric allele, OK = false 
					// 
					ad->ok[n] = ok;

					// 
					// save	SNPs to allele.
					// 
					//  	Note that first letter of non SNPs will also be saved but we
					//  	are guided by OK flag
					// 
					if ( ! ad->diallellic_snps[n] ) ad->diallellic_snps[n] = (char*)calloc(ncol+1, sizeof(char));
					for (k=0;k<ncol;k++)
					{
						ad->diallellic_snps[n][k] = buf[k];
					}

					// 
					// position along contig
					// 
					ad->loci[n] = Bp;
					n++;
				}
				nrow++;
			}
		}

		// number of diallelic sites
		ad->cnt_diallelic_snp_sites = n;
		printf ( "read %d sites / %d rows for %d cols from %s\n", ad->cnt_diallelic_snp_sites, nrow, ad->cnt_strains, filename );
		// maximum_locus position
		ad->maximum_locus = ad->loci[ad->cnt_diallelic_snp_sites-1];
		ad->map_pos_to_snp_index = (int*)calloc(ad->maximum_locus+1,sizeof(int));

		for (k=0;k<=ad->maximum_locus;k++)
			ad->map_pos_to_snp_index[k] = -1;

		// For non-missing all diallelic SNP sites
		// 		map_pos_to_snp_index[pos] = SNP index  
		// 			where pos = position along contig
		for (k=0;k<ad->cnt_diallelic_snp_sites;k++)
			if ( ad->ok[k] ) 
				ad->map_pos_to_snp_index[ad->loci[k]] = k;

		free(buffer);
		return(ad);
	}
}


ALLELE_DATA **ReadGenomeAlleleData( char *dir, char *maskfile, int nchr ) {

	int k, nseg=0;;
	char filename[256];

	// 
	// read allele data for each contig
	// 
	ALLELE_DATA** add = (ALLELE_DATA**)calloc(nchr+1,sizeof(ALLELE_DATA));
	// 
	// read GFF masking data for ignored regions
	// 
	SEGMENTS *mask = ReadGFF( maskfile, &nseg );
	for (k=1;k<=nchr;k++)
	{
		FILE *fp;
		sprintf(filename, "%s/chr%d.alleles.txt", dir, k);
		fp = fopen( filename, "r");
		if ( fp == NULL )
		{
			fprintf(stderr, "ERROR Could not open %s\n", filename );
			exit(1);
		}
		else
		{
			// chr   pse.bp   bp   nalleles   maf   bur-0   can-0   col-0   ct
			// 1     11       11   2          18    T       T       T       T 
			// 1     18       18   2          18    T       T       T       T 
			add[k] = ReadAlleleData( filename );
		}
	}

	// 
	//  set all masked regions to read as non OK
	// 
	// 
	for (k=0;k<nseg;k++)
	{
		int j;
		for (j=mask[k].from;j<=mask[k].to;j++)
		{
			add[mask[k].chr]->map_pos_to_snp_index[j] = -1;
		}
	}
	return(add);
}



//  Chr1	11	T	?
//	Chr1	18	T	?
//	Chr1	73	C	?
//	Chr1	83	T	?
SEQ_CALL *ReadSeqCall( char *filename, ALLELE_DATA **add, int nchr ) {
	FILE *fp = fopen(filename, "r");
	if ( NULL == fp )
	{
		fprintf(stderr, "ERROR Could not open %s\n", filename);
	}
	else
	{
		SEQ_CALL *sq = (SEQ_CALL*)calloc(1,sizeof(SEQ_CALL));
		int chr, bp, n=0, m=0;
		char r, a;

		// len = size of buffers
		sq->len=100000;
		sq->bp=(int*)calloc(sq->len,sizeof(int));
		sq->chr=(int*)calloc(sq->len,sizeof(int));
		sq->allele=(char*)calloc(sq->len,sizeof(char));
		sq->chr_start=(int*)calloc(1000,sizeof(int));
		sq->chr_end=(int*)calloc(1000,sizeof(int));
		sq->founder = (char**)calloc(sq->len,sizeof(char*));
		sq->mosaic = (SEGMENTS**)calloc(nchr+1,sizeof(SEGMENTS*));
		sq->nseg = (int*)calloc(nchr+1,sizeof(int));
		for (chr=1;chr<= nchr;chr++)
		{
			sq->chr_start[chr] = -1;
			sq->chr_end[chr] = -1;
		}

		// bp = pos
		// !!!!
		// assumes numeric chromosomes names
		// !!!!
		while (fscanf( fp, "Chr%d %d %c %c\n", &chr, &bp, &r, &a )==4 )
		{
			// ignore question marks
			// ignore if beyond maximum_locus known alleles
			if ( isalpha(a) && bp <= add[chr]->maximum_locus )
			{
				// 
				// look up founder reference allele from contig position
				// 		idx is offset into founder reference allele array
				//
				int idx = add[chr]->map_pos_to_snp_index[bp] ;
				m++;

				// 
				// ignore genotypes not in any of the founder strains
				// 
				if ( idx >= 0 )
				{
					if ( n >= sq->len )
					{
						sq->len *= 1.5;
						sq->bp=(int*)realloc(sq->bp, sq->len*sizeof(int));
						sq->chr=(int*)realloc(sq->chr, sq->len*sizeof(int));
						sq->allele=(char*)realloc(sq->allele, sq->len*sizeof(char));
						sq->founder = (char**)realloc(sq->founder, sq->len*sizeof(char*));
					}
					sq->bp[n] = bp;
					sq->chr[n] = chr;
					sq->allele[n] = a;
					// pointer array of alleles in founder strains corresponding to this contig position
					sq->founder[n] = add[chr]->allele[idx];

					// set contig bounds 
					if ( sq->chr_start[chr] == -1 || sq->chr_start[chr] > n ) sq->chr_start[chr] = n;
					if ( sq->chr_end[chr] < n )	sq->chr_end[chr] = n;

					// add founder strain names
					if ( n==0 )
					{
						sq->Nfounders = add[chr]->cnt_strains;
						sq->founder_name =add[chr]->strain_names;
					}

					n++;
				}
			}
		}

		fclose(fp);
		sq->len = n;
		//
		//  m = number of alleles within range /
		//  n = number of alleles within range and in founder strains
		// 
		printf( "red %d/%d alleles from %s\n", n, m, filename );
		return(sq);
	}
}

SEQ_CALL **ReadSeqCalls( char *dir, char *suffix, ALLELE_DATA **add, int nchr, int *nseq ) {
	glob_t globbuf;
	char pattern[256];
	struct dirent **d;
	sprintf(pattern, "%s/*%s", dir, suffix );
	glob( pattern, 0, NULL, &globbuf );
	int n = globbuf.gl_pathc;
	int chr;
	if (n < 0)
		perror("scandir");
	else
	{
		SEQ_CALL **seqs = (SEQ_CALL**)calloc(n,sizeof(SEQ_CALL));
		int k=0, penalty=100;
		while (k<n)
		{
			printf("%s\n", globbuf.gl_pathv[k]);
			seqs[k] = ReadSeqCall(  globbuf.gl_pathv[k], add, nchr ); 
			/*      for(chr=1;chr<=nchr;chr++) {
		  reconstruct_mosaic_diploid( seqs[k], penalty, chr );
		  } */
			k++;
		}
		*nseq = n;
		return(seqs);
	}
}


SEGMENTS *ReadGFF( char *filename, int *nsegments ) {

	FILE *fp  = fopen(filename, "r");

	if ( fp == NULL )
	{
		fprintf(stderr, "ERROR Could not open %s\n", filename);
	}
	else
	{
		int chr, from, to, nseg=100000, n=0;
		char build[256], type[256];
		SEGMENTS *segs = (SEGMENTS*)calloc(nseg,sizeof(SEGMENTS));
		while ( fscanf( fp, "Chr%d %s %s %d %d", &chr, build, type, &from, &to ) ==5 )
		{
			if ( n >= nseg )
			{
				nseg *= 1.5;
				segs = (SEGMENTS*)realloc( segs, sizeof(SEGMENTS)*nseg);
			}
			segs[n].chr = chr;
			segs[n].from = from;
			segs[n].to = to;
			n++;
		}
		*nsegments = n;
		return(segs);
	}
}

int reconstruct_mosaic( SEQ_CALL *sq, int penalty, int chr ) {

	int ns = sq->chr_end[chr]-sq->chr_start[chr]+1;
	int nf = sq->Nfounders;
	int offset = sq->chr_start[chr];
	int i, j, k, p, y;

	char **mat = (char**)calloc(ns,sizeof(char*));
	char **path = (char**)calloc(ns,sizeof(char*));
	int **opt = (int**)calloc(ns,sizeof(int*));
	int *step = (int*)calloc(ns+1,sizeof(int));
	int *jump = (int*)calloc(ns+1,sizeof(int));

	for (i=0;i<ns;i++)
	{
		mat[i] = (char*)calloc(nf,sizeof(char));
		path[i] = (char*)calloc(nf,sizeof(char));
		opt[i] = (int*)calloc(nf,sizeof(int));
	}

	for (i=0;i<ns;i++)
	{
		int k = i+offset;
		for (j=0;j<nf;j++)
			mat[i][j] = sq->allele[k] == sq->founder[k][j];
	}

	for (j=0;j<nf;j++)
	{
		opt[0][j] = mat[0][j];
		path[0][j] = j;
	}

	for (i=1;i<ns;i++)
	{
		for (j=0;j<nf;j++)
		{
			opt[i][j] = -1;
			path[i][j] = -1;
			for (k=0;k<nf;k++)
			{
				int x = opt[i-1][k] + mat[i][j] - penalty*(j!=k);
				if ( x > opt[i][j] )
				{
					opt[i][j] = x;
					path[i][j] = k;
				}
			}
		}
	}

	y = opt[ns-1][0];
	p = 0;
	for (j=0;j<nf;j++)
	{
		if ( opt[ns-1][j] > y )
		{
			y = opt[ns-1][j];
			p = j;
		}
	}

	step[ns] = p;
	i = ns;
	while (i>0)
	{
		p = path[--i][p];
		step[i] = p;
	} 

	jump[0] = 0;
	k = 1;
	for (i=1;i<=ns;i++)
	{
		if ( step[i] != step[i-1] )	jump[k++] = i;
	}
	jump[k] = ns;

	sq->mosaic[chr] = (SEGMENTS*)calloc(k,sizeof(SEGMENTS));
	for (i=0;i<k;i++)
	{
		int len = sq->bp[offset+jump[i+1]-1]-sq->bp[offset+jump[i]];
		int qlen = jump[i+1]-jump[i];
		double q, r;

		int m=0;

		for (j=jump[i];j<jump[i+1];j++)
			m += !mat[j][step[jump[i]]];
		q = qlen>0 ? m/(double)qlen : 0;
		r = len >0 ? m/(double)len : 0;
		j=jump[i];
		printf( "%d %s\t %10d %0d %10d %8d %8d %e %e\n", chr, sq->founder_name[step[j]], sq->bp[offset+j], sq->bp[offset+jump[i+1]-1], len, qlen, m, q, r );
		sq->mosaic[chr][k].chr = chr;
		sq->mosaic[chr][k].from = offset+j; 
		sq->mosaic[chr][k].to = offset+jump[i+1]-1;
		sq->mosaic[chr][k].id = step[j];
		sq->mosaic[chr][k].error = q;
	}
	sq->nseg[chr] = k;

	free(jump);
	free(step);

	for (i=0;i<ns;i++)
	{
		free(mat[i]);
		free(path[i]);
		free(opt[i]);
	}
	free(mat);
	free(path);
	free(opt);

	return(0);
}

int reconstruct_mosaic_diploid( SEQ_CALL *sq, int penalty, int chr ) {

	int ns = sq->chr_end[chr]-sq->chr_start[chr]+1;
	int nf = sq->Nfounders;
	int offset = sq->chr_start[chr];
	int g, h, i, j, k, p, y, nf2=nf*(nf+1)/2;
	char **name = (char**)calloc(nf2,sizeof(char*));

	char **mat = (char**)calloc(ns,sizeof(char*));
	short **path = (short**)calloc(ns,sizeof(short*));
	int **opt = (int**)calloc(ns,sizeof(int*));
	int *step = (int*)calloc(ns+1,sizeof(int));
	int *jump = (int*)calloc(ns+1,sizeof(int));

	for (i=0;i<ns;i++)
	{
		mat[i] = (char*)calloc(nf2,sizeof(char));
		path[i] = (short*)calloc(nf2,sizeof(short));
		opt[i] = (int*)calloc(nf2,sizeof(int));
	}

	for (i=0;i<ns;i++)
	{
		int k = i+offset;
		for (j=0,h=0;j<nf;j++)
		{
			for (g=0;g<=j;g++,h++)
			{
				mat[i][h] = ((sq->allele[k] == sq->founder[k][j]) || (sq->allele[k] == sq->founder[k][g]))*(1+ (j==g));
			}
		}
	}

	for (j=0,h=0;j<nf;j++)
	{
		for (g=0;g<=j;g++,h++)
		{
			name[h] = (char*)calloc(256, sizeof(char));
			sprintf(name[h], "%s %s %s", j==g? "=" : "*", sq->founder_name[j], sq->founder_name[g]);
			/*      printf("%d %d %d %s\n", j, g, h, name[h]);*/
		}
	}

	for (j=0;j<nf2;j++)
	{
		opt[0][j] = mat[0][j];
		path[0][j] = j;
	}

	for (i=1;i<ns;i++)
	{
		for (j=0;j<nf2;j++)
		{
			opt[i][j] = -1;
			path[i][j] = -1;
			for (k=0;k<nf2;k++)
			{
				int x = opt[i-1][k] + mat[i][j] - penalty*(j!=k);
				if ( x > opt[i][j] )
				{
					opt[i][j] = x;
					path[i][j] = k;
				}
			}
		}
	}

	y = opt[ns-1][0];
	p = 0;
	for (j=0;j<nf2;j++)
	{
		if ( opt[ns-1][j] > y )
		{
			y = opt[ns-1][j];
			p = j;
		}
	}

	step[ns] = p;
	i = ns;
	while (i>0)
	{
		p = path[--i][p];
		step[i] = p;
	} 

	jump[0] = 0;
	k = 1;
	for (i=1;i<=ns;i++)
	{
		if ( step[i] != step[i-1] )	jump[k++] = i;
	}
	jump[k] = ns;


	for (i=0;i<k;i++)
	{
		int len = sq->bp[offset+jump[i+1]-1]-sq->bp[offset+jump[i]];
		int qlen = jump[i+1]-jump[i];
		double q, r;
		int m=0;

		for (j=jump[i];j<jump[i+1];j++)
			m += !mat[j][step[jump[i]]];
		q = qlen>0 ? m/(double)qlen : 0;
		r = len >0 ? m/(double)len : 0;
		j=jump[i];

		printf( "%d %s\t %10d %0d %10d %8d %8d %e %e\n", chr, name[step[j]], sq->bp[offset+j], sq->bp[offset+jump[i+1]-1], len, qlen, m, q, r );
	}

	free(jump);
	free(step);

	for (i=0;i<ns;i++)
	{
		free(mat[i]);
		free(path[i]);
		free(opt[i]);
	}
	free(mat);
	free(path);
	free(opt);

	return(0);
}


